#setwd("filepath")
path<-getwd()
library(coda)
set.seed(666)

### IA court models and data
load("demsq_models.rda")
load("dif_models.rda")
load("iacrt_data_NAs.rda")

iacrt.models<-the.hi.model.list
iacrt.dif.models<-the.hi.model.dif.list

iacrt.dat<-newest.dat
dat<-na.omit(iacrt.dat)
dat$end<-ifelse(dat$year==2004,1,0)
dat$end<-ifelse(dat$iacrt==1,1,dat$end)
dat$id<-rev(cumsum(rev(dat$end)))
dat<-dat[order(dat$id, dat$year),]

### CCPR models and data
load("ccpr_demsq_models.rda")
load("ccpr_dif_models.rda")
load("ccpr_data_NAs.rda")

ccpr.models<-the.hi.model.list
ccpr.dif.models<-the.hi.model.dif.list

newest.dat<-newest.dat[!is.na(newest.dat$mean),]

### extract IA court model estimates

ests.list<-list()
for (i in 1:10){
	ests1<-as.matrix(iacrt.models[[i]][[1]])	
	ests2<-as.matrix(iacrt.models[[i]][[2]])	
	ests<-rbind(ests1, ests2)
	ests.list[[i]]<-ests
	}

ests.dif.list<-list()
for (i in 1:50){
	ests1<-as.matrix(iacrt.dif.models[[i]][[1]])	
	ests2<-as.matrix(iacrt.dif.models[[i]][[2]])	
	ests<-rbind(ests1, ests2)
	ests.dif.list[[i]]<-ests
	}

#### effect of democracy on IA court from democracy squared models
probs.list<-list()
dem<-seq(from=-1.44, to=2.01, by=.01)
for (i in 1:10){

	probs<-matrix(nrow=346, ncol=4)

	for (j in 1:346){
		pr<-
		plogis(
		ests.list[[i]][,34]+
		ests.list[[i]][,30]*dem[j]+
		ests.list[[i]][,32]*dem[j]^2+
		ests.list[[i]][,33]*20
		)
	
		probs[j,1]<-quantile(pr, probs=.05)
		probs[j,2]<-mean(pr)
		probs[j,3]<-quantile(pr, probs=.95)
		probs[j,4]<-dem[i]
		}

	probs.list[[i]]<-probs

	}
	
### extract CCPR model estimates

ests.list.2<-list()
for (i in 1:10){
	ests1<-as.matrix(ccpr.models[[i]][[1]])	
	ests2<-as.matrix(ccpr.models[[i]][[2]])	
	ests<-rbind(ests1, ests2)
	ests.list.2[[i]]<-ests
	}

ests.dif.list.2<-list()
for (i in 1:50){
	ests1<-as.matrix(ccpr.dif.models[[i]][[1]])	
	ests2<-as.matrix(ccpr.dif.models[[i]][[2]])	
	ests<-rbind(ests1, ests2)
	ests.dif.list.2[[i]]<-ests
	}

#### substantive effect of democracy on CCPR from democracy squared models
probs.list.2<-list()
dem<-seq(from=-1.44, to=2.01, by=.01)
for (i in 1:10){

	probs<-matrix(nrow=346, ncol=4)

	for (j in 1:346){
		pr<-
		plogis(
		ests.list.2[[i]][,33]+
		ests.list.2[[i]][,29]*dem[j]+
		ests.list.2[[i]][,31]*dem[j]^2+
		ests.list.2[[i]][,32]*20
		)
	
		probs[j,1]<-quantile(pr, probs=.05)
		probs[j,2]<-mean(pr)
		probs[j,3]<-quantile(pr, probs=.95)
		probs[j,4]<-dem[i]
		}

	probs.list.2[[i]]<-probs

	}

tr.gr<-rgb(150,150,150,150,maxColorValue=255)

pdf("iacrt_demsq.pdf", width=11, height=5)

#dev.new(width=11, height=5)
par(mfrow=c(1,2), mar=c(4, 4, 2, 1), cex.lab=1.25,cex.main=1.5)
plot(density(iacrt.dat$mean), ylim=c(0, 1), xlim=c(-2, 2.5), main="Inter-American Court", xlab="Latent Democracy", ylab="Pr(Acceptance)", axes=F, type="n")
abline(h=c(0,0.2,0.4,0.6,0.8,1), col=tr.gr, lty=3, lwd=1.5)
polygon(density(iacrt.dat$mean), col="grey")
axis(1, at=c(-2, -1, 0, 1, 2))
axis(2, at=c(0, 0.2, 0.4, 0.6, 0.8, 1))
polygon(density(iacrt.dat$mean), col="grey")
for (i in 1:10){
	lines(dem, probs.list[[i]][,2], col="black", lty=1, lwd=3)
	lines(supsmu(dem, probs.list[[i]][,1]), col="gray48", lty=1, lwd=3)
	lines(supsmu(dem, probs.list[[i]][,3]), col="gray48", lty=1, lwd=3)
	}
rug(iacrt.dat$mean[iacrt.dat$iacrt==1], lwd=2)
box()

plot(density(newest.dat$mean), ylim=c(0, 1), xlim=c(-2, 2.5), main="ICCPR", xlab="Latent Democracy", ylab="Pr(Ratification)", axes=F, type="n")
abline(h=c(0,0.2,0.4,0.6,0.8,1), col=tr.gr, lty=3, lwd=1.5)
polygon(density(newest.dat$mean), col="grey")
axis(1, at=c(-2, -1, 0, 1, 2))
axis(2, at=c(0, 0.2, 0.4, 0.6, 0.8, 1))
polygon(density(newest.dat$mean), col="grey")
for (i in 1:10){
	lines(dem, probs.list.2[[i]][,2], col="black", lty=1, lwd=3)
	lines(supsmu(dem, probs.list.2[[i]][,1]), col="gray48", lty=1, lwd=3)
	lines(supsmu(dem, probs.list.2[[i]][,3]), col="gray48", lty=1, lwd=3)
	}
rug(newest.dat$mean[newest.dat$iccpr_rat==1], lwd=2)
box()

dev.off()

#### effect of transition for IA court (not calculated for CCPR)

probs.list<-list()
for (i in 1:50){

	probs<-matrix(nrow=30, ncol=4)

	for (j in 1:30){
		pr1<-
		plogis(
		ests.dif.list[[i]][,34]+
		ests.dif.list[[i]][,30]*0+
		ests.dif.list[[i]][,31]*0+
		ests.dif.list[[i]][,32]*j
		)
		
		pr2<-
		plogis(
		ests.dif.list[[i]][,34]+
		ests.dif.list[[i]][,30]*1+
		ests.dif.list[[i]][,31]*0+
		ests.dif.list[[i]][,32]*j
		)
		
		pr<-pr2-pr1
		probs[j,1]<-quantile(pr, probs=.05)
		probs[j,2]<-mean(pr)
		probs[j,3]<-quantile(pr, probs=.95)
		probs[j,4]<-j

		}

	probs.list[[i]]<-probs

	}

lab1<-expression(Democracy^2)
lab2<-expression(Delta[t-1])
lab3<-expression(Delta[t-2])
lab4<-expression(Delta[t-3])
lab5<-expression(Delta[t-4])
lab6<-expression(Delta[t-5])
dif.labs<-c(lab2,lab3,lab4,lab5,lab6)

probs.t25<-matrix(nrow=length(probs.list),ncol=4)
for (i in 1:length(probs.list)){
	probs.t25[i,]<-probs.list[[i]][20,]
	}

pdf("iacrt_dif_20.pdf")
par(mfrow=c(5,1), mar=c(2, 2, 2, 1))
for (j in 1:5){
	plot(seq(1:10), probs.t25[((j*10)-(9)):(j*10),2], type="p", pch=16, cex=1.75, ylab="", xlab="", ylim=c(-0.05,0.6), main=dif.labs[j], cex.main=1.5, axes=F)
	for (i in 0:9){
		lines(c(i+1, i+1), c(probs.t25[((j*10)-(9))+i,1],probs.t25[((j*10)-(9))+i,3]), lwd=1.75)
		}
	axis(1, at=seq(1:10), labels=c("","","","","","","","","",""))
	axis(2, at=c(0, 0.2, 0.4, 0.6), labels=c("0", "0.2", "0.4", "0.6"))
	abline(h=0, lty=2)
	abline(h=c(0.2,0.4,0.6), lty=3, col=tr.gr, lwd=1.5)
	box()
	}
dev.off()	

#### Time series plots
cc<-c(90,91,95,155,160,165)
pdf("iacrt_ts1.pdf")
par(mfrow=c(3, 2), mar=c(3, 2, 4, 1))
for (i in 1:length(cc)) {
#	tmp<-paste(as.character(cc[i]),"pdf",sep=".")
#	pdf(tmp)
	plot(iacrt.dat$mean[iacrt.dat$ccode==cc[i]]~iacrt.dat$year[iacrt.dat$ccode==cc[i]], type="l", lwd=2, xlab="",ylab="", main=iacrt.dat$country.x[iacrt.dat$ccode==cc[i]][1], ylim=c(-1.75, 2.5))	
	polygon(c(iacrt.dat$year[iacrt.dat$ccode==cc[i]],rev(iacrt.dat$year[iacrt.dat$ccode==cc[i]]),iacrt.dat$year[iacrt.dat$ccode==cc[i]][1]),c(iacrt.dat$pct975[iacrt.dat$ccode==cc[i]],rev(iacrt.dat$pct025[iacrt.dat$ccode==cc[i]]),iacrt.dat$pct025[iacrt.dat$ccode==cc[i]][1]),col=tr.gr,border=NA)
	abline(v=iacrt.dat$year[iacrt.dat$ccode==cc[i]&iacrt.dat$iacrt==1][1], lty=2)
	abline(h=c(-1,0,1,2), lty=3, col=tr.gr, lwd=1.5)
#	dev.off()
}
dev.off()

iacrt.dat$country.y<-as.character(iacrt.dat$country.y)
iacrt.dat$country.y[iacrt.dat$cc==2]<-"USA"
iacrt.dat$country.y[iacrt.dat$cc==42]<-"DR"
cc<-setdiff(unique(dat$ccode),cc)
cc2<-cc[1:9]
cc3<-cc[10:23]
cc<-cc2
pdf("iacrt_ts2.pdf")
par(mfrow=c(3, 3), mar=c(3, 2, 3, 1))
for (i in 1:length(cc)) {
	plot(iacrt.dat$mean[iacrt.dat$ccode==cc[i]]~iacrt.dat$year[iacrt.dat$ccode==cc[i]], type="l", lwd=2, xlab="",ylab="", main=iacrt.dat$country.y[iacrt.dat$ccode==cc[i]][1], ylim=c(-1.75, 2.5))	
	polygon(c(iacrt.dat$year[iacrt.dat$ccode==cc[i]],rev(iacrt.dat$year[iacrt.dat$ccode==cc[i]]),iacrt.dat$year[iacrt.dat$ccode==cc[i]][1]),c(iacrt.dat$pct975[iacrt.dat$ccode==cc[i]],rev(iacrt.dat$pct025[iacrt.dat$ccode==cc[i]]),iacrt.dat$pct025[iacrt.dat$ccode==cc[i]][1]),col=tr.gr,border=NA)
	abline(v=iacrt.dat$year[iacrt.dat$ccode==cc[i]&iacrt.dat$iacrt==1][1], lty=2)
	abline(h=c(-1,0,1,2), lty=3, col=tr.gr, lwd=1.5)
}
dev.off()

cc<-cc3
pdf("iacrt_ts3.pdf")
par(mfrow=c(4, 4), mar=c(3, 2, 3, 1))
for (i in 1:length(cc)) {
	plot(iacrt.dat$mean[iacrt.dat$ccode==cc[i]]~iacrt.dat$year[iacrt.dat$ccode==cc[i]], type="l", lwd=2, xlab="",
	ylab="", main=iacrt.dat$country.y[iacrt.dat$ccode==cc[i]][1], ylim=c(-1.75, 2.5))	
	polygon(c(iacrt.dat$year[iacrt.dat$ccode==cc[i]],rev(iacrt.dat$year[iacrt.dat$ccode==cc[i]]),iacrt.dat$year[iacrt.dat$ccode==cc[i]][1]),c(iacrt.dat$pct975[iacrt.dat$ccode==cc[i]],rev(iacrt.dat$pct025[iacrt.dat$ccode==cc[i]]),iacrt.dat$pct025[iacrt.dat$ccode==cc[i]][1]),col=tr.gr,border=NA)
	abline(v=iacrt.dat$year[iacrt.dat$ccode==cc[i]&iacrt.dat$iacrt==1][1], lty=2)
	abline(h=c(-1,0,1,2), lty=3, col=tr.gr, lwd=1.5)
}
dev.off()

#### Ropeladder plots

### coefficient estimates from IA court models
dem.ests<-matrix(nrow=10,ncol=3)
demsq.ests<-matrix(nrow=10,ncol=3)
for (i in 1:length(ests.list)){
	dem.ests[i,1]<-mean(ests.list[[i]][,30])
	dem.ests[i,2:3]<-quantile(ests.list[[i]][,30],c(0.05,0.95))
	demsq.ests[i,1]<-mean(ests.list[[i]][,32])
	demsq.ests[i,2:3]<-quantile(ests.list[[i]][,32],c(0.05,0.95))
	}

dif.ests<-matrix(nrow=50,ncol=3)
for (i in 1:length(ests.dif.list)){
	dif.ests[i,1]<-mean(ests.dif.list[[i]][,30])
	dif.ests[i,2:3]<-quantile(ests.dif.list[[i]][,30],c(0.05,0.95))
	}

coefs<-matrix(nrow=9,ncol=3)
coefs[1,1]<-mean(dif.ests[41:50,1])
coefs[1,2]<-min(dif.ests[41:50,2])
coefs[1,3]<-max(dif.ests[41:50,3])

coefs[2,1]<-mean(dif.ests[31:40,1])
coefs[2,2]<-min(dif.ests[31:40,2])
coefs[2,3]<-max(dif.ests[31:40,3])

coefs[3,1]<-mean(dif.ests[21:30,1])
coefs[3,2]<-min(dif.ests[21:30,2])
coefs[3,3]<-max(dif.ests[21:30,3])

coefs[4,1]<-mean(dif.ests[11:20,1])
coefs[4,2]<-min(dif.ests[11:20,2])
coefs[4,3]<-max(dif.ests[11:20,3])

coefs[5,1]<-mean(dif.ests[1:10,1])
coefs[5,2]<-min(dif.ests[1:10,2])
coefs[5,3]<-max(dif.ests[1:10,3])

coefs[7,1]<-mean(demsq.ests[,1])
coefs[7,2]<-min(demsq.ests[,2])
coefs[7,3]<-max(demsq.ests[,3])
coefs[8,1]<-mean(dem.ests[,1])
coefs[8,2]<-min(dem.ests[,2])
coefs[8,3]<-max(dem.ests[,3])

### coefficient estimates for democracy squared models
dem.ests.2<-matrix(nrow=10,ncol=3)
demsq.ests.2<-matrix(nrow=10,ncol=3)
for (i in 1:length(ests.list.2)){
	dem.ests.2[i,1]<-mean(ests.list.2[[i]][,29])
	dem.ests.2[i,2:3]<-quantile(ests.list.2[[i]][,29],c(0.05,0.95))
	demsq.ests.2[i,1]<-mean(ests.list.2[[i]][,31])
	demsq.ests.2[i,2:3]<-quantile(ests.list.2[[i]][,31],c(0.05,0.95))
	}

dif.ests.2<-matrix(nrow=50,ncol=3)
for (i in 1:length(ests.dif.list.2)){
	dif.ests.2[i,1]<-mean(ests.dif.list.2[[i]][,29])
	dif.ests.2[i,2:3]<-quantile(ests.dif.list.2[[i]][,29],c(0.05,0.95))
	}

coefs.2<-matrix(nrow=9,ncol=3)
coefs.2[1,1]<-mean(dif.ests.2[41:50,1])
coefs.2[1,2]<-min(dif.ests.2[41:50,2])
coefs.2[1,3]<-max(dif.ests.2[41:50,3])

coefs.2[2,1]<-mean(dif.ests.2[31:40,1])
coefs.2[2,2]<-min(dif.ests.2[31:40,2])
coefs.2[2,3]<-max(dif.ests.2[31:40,3])

coefs.2[3,1]<-mean(dif.ests.2[21:30,1])
coefs.2[3,2]<-min(dif.ests.2[21:30,2])
coefs.2[3,3]<-max(dif.ests.2[21:30,3])

coefs.2[4,1]<-mean(dif.ests.2[11:20,1])
coefs.2[4,2]<-min(dif.ests.2[11:20,2])
coefs.2[4,3]<-max(dif.ests.2[11:20,3])

coefs.2[5,1]<-mean(dif.ests.2[1:10,1])
coefs.2[5,2]<-min(dif.ests.2[1:10,2])
coefs.2[5,3]<-max(dif.ests.2[1:10,3])

coefs.2[7,1]<-mean(demsq.ests.2[,1])
coefs.2[7,2]<-min(demsq.ests.2[,2])
coefs.2[7,3]<-max(demsq.ests.2[,3])
coefs.2[8,1]<-mean(dem.ests.2[,1])
coefs.2[8,2]<-min(dem.ests.2[,2])
coefs.2[8,3]<-max(dem.ests.2[,3])

ylabels1<-c(lab6, lab5, lab4, lab3, lab2)
lab0<-expression(Democracy)
ylabels2<-c(lab1, lab0)
ylabels4<-c("Models 2-6 ", " Model 1 ")

pdf(file="iacrt_coefs.pdf",)

par(mfrow=c(1,2), oma=c(1, 4, 1, 1), mar=c(2, 1.5, 1.5, 1), cex.axis=0.85, cex.main=1.25)
plot(coefs[,1], seq(1:dim(coefs)[1]), type="p", pch=16, cex=1.5, ylab="", xlab="", axes=F, xlim=c(-5, 7.5), main="Inter-American Court")
for (i in 1:dim(coefs)[1]){
	lines(coefs[i,2:3], c(i, i), lwd=1.5)
}
axis(2, labels=ylabels1, at=seq(1:5), las=1)
axis(2, labels=ylabels2, at=c(7, 8), las=1)
axis(2, labels=ylabels4, at=c(6, 9), las=1, font=2, tick=F)
axis(1, at=seq(from=-4, to=6, by=2))
abline(v=0, lty=2)
abline(v=c(-4,-2,2,4,6), lty=3, col=tr.gr, lwd=1.5)
box()

plot(coefs.2[,1], seq(1:dim(coefs.2)[1]), type="p", pch=16, cex=1.5, ylab="", xlab="", axes=F, xlim=c(-5, 7.5), main="ICCPR")
for (i in 1:dim(coefs.2)[1]){
	lines(coefs.2[i,2:3], c(i, i), lwd=1.5)
}
axis(1, at=seq(from=-4, to=6, by=2))
abline(v=0, lty=2)
abline(v=c(-4,-2,2,4,6), lty=3, col=tr.gr, lwd=1.5)
box()

dev.off()

